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ABSTRACT 



> 

Q\ ; Planetary embryos embedded in a gas disc suffer a decay in semimajor axis - 

type I migration - due to the asymmetric torques produced by the interior and 
exterior wakes raised by the body (Goldreich & Tremaine 1980; Ward 1986). This 
' presents a challenge for standard oligarchic approaches to forming the terrestrial 

planets (Kokubo & Ida 1998) as the timescale to grow the progenitor objects near 
1 AU is longer than that for them to decay into the Sun. In this paper we inves- 
tigate the middle and late stages of oligarchic growth using both semi-analytic 
methods (based upon Thommes et al. 2003) and N-body integrations, and vary 
gas properties such as dissipation timescale in different models of the protoplan- 
etary disc. We conclude that even for near-nominal migration efficiencies and 
gas dissipation timescales of ~ 1 Myr it is possible to maintain sufficient mass in 
?_i ■ the terrestrial region to form Earth and Venus if the disc mass is enhanced by 



factors of ~ 2 — 4 over the minimum mass model. The resulting configurations 
differ in several ways from the initial conditions used in previous simulations of 
the final stages of terrestrial accretion (e.g. Chambers 2001), chiefly in (1) larger 
inter-embryo spacings, (2) larger embryo masses, and (3) up to ~0.4M e of ma- 
terial left in the form of planetesimals when the gas vanishes. The systems we 
produce are reasonably stable for ~ 100 Myr and therefore require an external 
source to stir up the embryos sufficiently to produce final systems resembling the 
terrestrial planets. 



Subject headings: planets, formation 
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1. Introduction 

A standard model for forming the terrestrial planets has emerged, and is divided broadly 
into three stages (see references in Canup et al. 2000). The first stage is poorly understood, 
but the formation of the Sun is thought to leave a large but thin protoplanetary disc of 
both gas and solids, with coagulation or gravitational instability producing rocky objects 
('planetesimals') in the metre-to-km size range. In the second stage, interactions between 
the planetesimals result in an early phase of 'runaway growth', to be discussed below, where 
the single largest body in a given region breaks away from the mass spectrum and becomes 
significantly larger than the second-largest. Over ~ 10 5 — 10 6 years, this gives rise to an 
oligarchic phase, also discussed below, which produces many well-spaced comparable-mass 
embryos. In the third stage, the embryos interact with each other and what remains of the 
the disc, merging and producing the terrestrial planets we observe today on a timescale of 
~ 10 7 - 10 8 years. 

The basic paradigm for the second stage, the 'planetesimal problem', is described in 
Wetherill & Stewart (1993). One considers a very large number of small bodies of compara- 
ble mass, moving in low-eccentricity, low-inclination orbits around the Sun. These objects 
will suffer close encounters with each other, changing the velocity distribution, and occa- 
sionally collide and either merge or fragment, changing the mass spectrum as well. Energy 
equipartition in the system decreases the random velocities of the larger bodies, and increases 
the velocities of the smaller ones, an effect known as 'dynamical friction'. It has been shown 
(Wetherill & Stewart 1989) that dynamical friction leads to runaway growth of the largest 
body, as the decrease in velocity increases the collisional cross section and thus the accretion 
rate. Later simulations (e.g. Ida & Makino 1992a,b; Kokubo & Ida 1996, 1998) confirmed 
this, and demonstrated that the end result of runaway growth was an 'oligarchic' phase where 
the accretion process is dominated by the gravitational stirring of well-separated embryos of 
roughly equal masses. 

All simulations published to date of the second stage have either neglected any interac- 
tion between the gas disk and the solid bodies or have included only the aerodynamic gas 
drag acting on the solid bodies as described by Adachi et al. (1976). However, it is well 
known that the interactions between the protoplanetary disc and the protoplanets can also 
cause the latter to migrate. This can happen in several ways: type I, where the migration is 
driven by the asymmetry between the torques generated by the interior and exterior wakes 
of the body (Goldreich & Tremaine 1980; Ward 1986); type II, where the protoplanet is large 
enough to open up a gap in the disc, coupling it to the viscous evolution of the disc (Pa- 
paloizou & Lin 1984); and recently type III (Masset & Papaloizou 2003; Artymowicz 2004), 
a very fast migration mode which occurs when the embryo begins to migrate quickly enough 
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that corotating material cannot librate. (The solid material in the disc can also be scattered 
by the embryos, leading to embryo migration directly; this 'type 0' migration (Fernandez & 
Ip 1984) is thought to play an important role in the evolution of the outer solar system.) 

In most models, protoplanets of terrestrial mass are incapable of opening a gap and 
therefore do not suffer type II migration, and likewise type III should not be relevant (Masset 
& Papaloizou 2003). Accordingly we restrict ourselves to considering the effects of type 
I migration, specifically the tidal interactions between the planet and disc which result in 
semimajor axis decay and random velocity decay (both eccentricity and inclination damping.) 
In the models we are considering, the typical effect is to decrease the semimajor axis, which 
leads to difficulties for accretion: the type I migration rate scales with the mass of the object, 
so as an object grows it moves inwards faster. Indeed, the type I migration timescales for 
the embryos produced near 1 AU are typically comparable to the timescales for their growth; 
this problem is discussed in Tanaka et al. (2002). 

In a sense, we seek to refine the initial conditions for the third stage, by including the 
effects of the tidal interaction during the second stage. This is very similar in spirit to 
Kominami & Ida (2004), but they include only the effects of random velocity damping, not 
of semimajor axis decay, and study only the third stage itself. To this end we will investigate 
the effects of different assumptions about the disc and the migration efficiency on the state of 
the protoplanetary system when the gas vanishes, and thus our ability to produce terrestrial 
planets from what material remains. Note that we concentrate here on forming the Earth 
and Venus. Mercury has a small enough mass (0.055 M e ) that it is easily considered debris 
at this scale, and the problem of forming Mars has its own special difficulties involving the 
asteroid belt, the possible presence of other large embryos, and influences from Jupiter, which 
we defer. We will not discuss here the third stage resulting from this work, because it is 
probably dependent on the details of the stirring of the embryos near the end of the second 
stage. This is likely to involve the perturbations introduced by the formation of Jupiter 
and Saturn, for which models are sufficiently uncertain that we reserve a discussion of these 
effects for a subsequent paper. 

In §2 we develop a semi-analytic model for oligarchy during type I migration in the ter- 
restrial region and use the results to motivate the N-body simulations which follow. In §3 we 
discuss our numerical N-body methods, and present the results of our N-body simulations in 
§4. Section 5 contains a further discussion of our scenario and we summarize our conclusions 
in §6. 
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2. Models 

In this section we construct a model of oligarchic formation by adapting the model of 
Thommes et al. (2003), which built on the earlier work of Kokubo & Ida (1998, 2000, 2002). 
We content ourselves with summarizing their results and noting our modifications. We then 
discuss the differences between the physics in this semi-analytic model and the physics we 
use in our N-body simulations. 

We use the semi-analytic model for two distinct purposes: first, to provide a quick 
way to explore parameter space and determine what regimes are good candidates for more 
detailed N-body investigation; and second, to generate initial conditions for said N-body 
simulations. This allows us to begin our runs later in the process than would otherwise be 
possible, thereby avoiding the expense of handling very large numbers of embryos at early 
times. 



2.1. Semi-analytic Model 

Thommes et al. (2003) represent the oligarchic system by a two-component model con- 
sisting of embryos and field planetesimals, where the former accrete but the latter do not. 
The model variables are the embryo mass M and planetesimal surface density S m as func- 
tions of semimajor axis a, and we seek to evolve M(a) and S m (a) over time t. 

Note that since M(a) is continuous but is meant to represent a discrete population, 
it corresponds not to the mass of material in the form of embryos in a given region, but 
to the mass that an embryo would have were one present at that semimajor axis. It is 
also assumed that the embryo separation b measured in single-planet Hill units (defined by 
Th = (M/3M ) 1 / 3 a) is constant and contains all the information about the effects of any 
embryo-embryo mergers in the model. This assumption is justified in Kokubo & Ida (1998) 
as the consequence of an equilibration between the decrease in b caused by accretion of small 
objects and an increase in b caused by two-body scattering followed by eccentricity damping 
due to dynamical friction from the field of planetesimals. The properties of the resulting 
(non-continuous) embryo population can be determined from M (a) by discretizing the curve 
appropriately. Thus we construct an effective embryo surface density Em = M/(2naAa), 
setting the embryo spacing Aa = bru- 

The embryos are assumed to have uniform physical density pu- It is also assumed that 
the embryos are kept in near-circular orbits by dynamical friction with the field (or, when 
active, tidal damping from the gas) and therefore the important random velocity is that of 
the field particles. 
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The field population is similarly approximated as being comprised of planetesimals of 
uniform mass m and density p m (and therefore uniform radius r m ) with (root-mean-square) 
eccentricity e m and inclination i m satisfying e m = 2i m . Thommes et al. (2003) determine 
the eccentricity, after Ida & Makino (1993), by equating the stirring timescale of the field 
by the embryos with the damping timescale due to aerodynamic drag on the field particles. 

The stirring timescale in a dispersion-dominated regime can be written as 



stir — 



M 



40 V CM / £ M a 2 lT 



where G is the gravitational constant and Q the orbital frequency. We use the damping 
timescale due to aerodynamic drag from Adachi et al. (1976), where p gas is the density of 
the gas and Co is a drag efficiency factor (~ 1) dependent upon the Reynolds number of the 
disc. We neglect resonant interactions. In our notation, 

1 Tfl 

aero = 4^/2)^2^^ ( } 

Equating these two timescales yields an equilibrium eccentricity of 

1.72 m 1/15 M 1/3 pm /15 

Thommes et al. (2003) derive an expression for the growth of the embryo mass 

dM _ 3.93 il4 /6 G l/2 S m M 2 / 3 Cf p 2 £ 



dt p V3 a i/io m 2/i 5 ^/i5 



(4) 



and a corresponding decrease in planetesimal surface density from conservation of mass 

dX m _ -M'J 3 dM 
dt 3 2 / 3 67ra 2 MV3 dt ■ U 

(The aerodynamic drag also affects E m , as we discuss later.) 

We restrict consideration of initial disc conditions at T = to those resembling the 
minimum mass model (Hayashi 1981). Specifically, we assume the surface density in solids 
has the form S solid = /c(r/AU) p ; the nominal values are k = 7.1gcm~ 2 ,p = —1.5 and 
we introduce an enhancement factor f cnh so that S so i id = / enh (7.1 gem -2 ) (r/AU) p . We 
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adopt a gas-to-solid ratio of 240, so that £ gas = 240 E S oiid- We assume the disc has an 
exponential vertical structure such that p gas (r, z) = p gas ( r ) exp(— z 2 /zq), with p gas ( r ) taken 
at the midplane and Zq(t) the disc half-thickness. (Integrating the z-dependent term over 
all z gives the relation p gas = E gas / \Znz .) The half-thickness varies as z (r) = ^(r/AU) 5 / 4 , 
where we take the thickness at 1 AU to be Z\ = 0.07 AU. We model the dissipation of the gas 
by introducing an exponential decay with a characteristic time Td ccay which is fixed for a given 
model and does not vary through the disc, p oc exp(— t/r decay ), where we allow Td ecay = oo. 
We define a parameter r\ which measures the degree of pressure support of the disc (and 
thus deviations of circular orbits from Keplerian velocity): rj = (n/16)(a + (3)(z /a) 2 where 
a = 5/4 — p and (3 = 1/2; typically rj ~ 0.001. We assume the gas is in cylindrical rotation 
such that t> gas = fkcpVl ~~ 277. 

For the planetesimal migration introduced by aerodynamic drag, we adopt the orbit- 
averaged approximation valid for small e, i, and rj due to Adachi et al. (1976): 



= ft L * ~ 2 T±r m (i* + + " 2 ) " 2 {" + (7 + te) < + (6) 

We also add a prescription for type I migration based on that of Papaloizou & Larwood 
(2000), assuming that the embryos are on circular orbits (see also Tanaka et al. 2002). They 
derive a timescale for semimajor axis damping t a : 



1 / a 3 \ 1/2 (Z±_ /_a_x i/A 2 / S gas vra 2 \ ( M_\ 1 
ta c a \GM Q J Uu VAUv 1 M M Q J \m g ) ' Uj 
where c a is a migration efficiency parameter, with c a — 1 being nominal and c a = yielding 
t a = 00. By varying c a in our simulations we can account for an uncertainty of order unity 
in the rate of migration. 

For the purpose of this model we can construct a rough radial velocity for protoplanet 
material 



da 

VM = Tt 



= -- (8) 

tidal t n 



We emphasize that since it is not clear (to choose one assumption among many) whether 
the aerodynamic drag coefficient C D should be ~1 or ~2 (Adachi et al. 1976), the above 
model is only expected to be accurate to within factors of order unity. This is true even if 
we neglect consequences of assuming that all field particles have the same fixed mass and 
the issue of the evolution of inter-embryo spacing (see §4.4). 
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2.2. Results from Semi-analytic Model 

We now integrate the model described in the previous section forward in time. We have 
two variables, the embryo mass M(a) and the planetesimal surface density S m (a), whose 
evolution (in the absence of aerodynamic or tidal migration) is specified by equations (4) 
and (5). We introduce migration through equations (6), (7), and (8). These five equations 
are simultaneously integrated across a zone from a = 0.2 AU to a = 5.0 AU, subject to our 
time-dependent gas profile. 

Figure 1 shows the evolution of the standard minimum-mass model (k = 7.1 gem" 2 , 
/enh = 1-0, p = -1.5, r dccay = oo) assuming that b = 10. We arbitrarily set M = 10 _6 M e 
as the seed mass for the embryos at all a, so that the total fraction of mass in embryos is 
< 5 x 10~ 4 that of the mass in the planetesimal field; accretion at early times is so quick 
that the precise value is unimportant as long as M <C Mfi na i. As is evident on the left part 
of the curves, the embryo mass nearly asymptotes toward a power-law in semimajor axis a: 
except for deviations introduced by planetesimal migration and the nonzero initial mass M 
of the embryo, the limiting mass is set solely by the amount of planetesimal material and 
the chosen b, M\ im = A/8/(3M )(S m nba 2 ) 3 ^ 2 . (The right side of the curve will reach the 
same asymptote eventually.) 

Figure 2 shows the minimum-mass model with tidal migration added (keeping Td eC ay = 
oo): several differences from the previous case are immediately apparent. The maximum 
embryo mass is substantially reduced, from 0.086M e at 10 Myr to 0.027M e . Without 
migration there is a sharp drop-off in embryo mass past the peak of the accretion front, but 
with tidal migration there is a wide band in semimajor axis where the embryo mass is within 
25% of the maximum. With Td eC a y = oo there is nothing to prevent embryo material from 
continuing to migrate inwards as it accretes, and thus eventually all embryo mass falls into 
the Sun. 

Figure 3 is the result of introducing gas dissipation of r decay = 1 Myr, which produces an 
intermediate system. The maximum embryo mass reached is ~0.043M e , down by a factor 
of 2 from the no-migration case, but the gas decay has several effects: type I migration is 
halted, preserving embryo mass; the increase in planetesimal eccentricity as a result of the 
decrease in aerodynamic drag increases the accretion time; and as a result of this increase, 
a fair amount of material is left in the form of planetesimals. However, the total amount of 
material in the terrestrial region (0.5AU < a < 1.5AU) is insufficient to later accrete into 
the Earth and Venus during stage three. 

The dissipation timescale of ~ 1 Myr was chosen based on studies of the disc lifetimes 
around stars in young clusters by Haisch et al. 2001. They find that roughly half the stars 
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lose their discs in < 3 Myr, and the overall disc lifetime was ~ 6 Myr. (In the N-body 
simulations to be discussed later, we will consider discs with both larger and greater Td eca y) 

This suggests that in order to obtain the appropriate amount of material in the terrestrial 
region with nominal migration, one could (1) introduce a dissipating disc and (2) enhance 
the disc above the minimum mass model. Figure 4 compares the end state after 100 Myr 
of evolution for the minimum mass model with various Td eC a y values with the same for a 
disc enhanced by a factor of three. The enhanced case both moves the point of maximum 
embryo mass outwards and raises the maximum mass reached by almost a factor of ~ 3, but 
is otherwise qualitatively similar to the minimum mass results. From this and other runs 
(not shown) we conclude that if we believe the nominal migration rates, then to within the 
expected reliability of the model, disc enhancements of a factor of several and dissipation 
times of ~1 Myr should suffice to keep enough prototerrestrial mass in the region. 

Unfortunately these analytic models are sufficiently simplified that they miss many in- 
teresting dynamical interactions and so cannot predict the detailed behaviour of the embryos. 
This is especially true at the later stages when the embryo eccentricities are no longer being 
damped by the gas. The assumption that 6~10 also breaks down (see §4.4) as the embryo 
mass increases, and therefore semi-analytic estimates of the mass in embryos remaining at 
Myr timescales may be unreliable. We therefore want to construct N-body realizations of 
these models and integrate them directly, as we discuss next. 



2.3. N-body modifications to Semi-analytic model 

For consistency, we carry as much of the above approach as possible through to our 
N-body code. There are nevertheless modifications in both the aerodynamic and the tidal 
drag formulae. 

First, instead of the approximation for the aerodynamic drag in eq. (6), we consider the 
planetesimal migration induced by aerodynamic drag 



_ dv _ V ~ Vgas , v 

"aero j \J ) 

7~acro 

where v is the Cartesian velocity of the object, v gas the velocity of the gas at the object's 
position, and r acro is given by eq. (2). As before, we assume that v gas = v^y/l ~~ ^V- 

Second, for the tidal drag, we use the full approach of Papaloizou & Larwood (2000), 
which the authors developed to handle the case where a protoplanet's eccentricity can be 



minimum mass model at T = 1 00 Myr 




x decay = 0.1 Myr 

tdecay = 1 Myr 
^decay = 1 Myr 



semimajor axis [AU] 



greater than the scale height-to-semimajor axis ratio. They derive timescales for semimajor 
axis damping t a and for eccentricity damping t e : 



c a \GMcJ \A\J \A\JJ J \ M J 
'My 1 / l + (e/1.3) (VAIJ)' 5 (a/AU)' 5/4 \ 
Mj V 1 - (e/l-l) (^/AU)-4 ( a /AU)-i J ' 



e c e \GM Q J \axj \axj) ) \ m & 



1 + 7 aTF TTT h (11) 



M 7 \ 4 \AUy VAU7 

where (as before) c a is the migration efficiency, and c e is the analogous damping efficiency. 

They also argue that if the inclination damping timescale is not significantly shorter 
than the eccentricity damping timescale then it plays little role in the equilibrium state; we 
set ti = t e for simplicity. From these we can find the acceleration on an object due to tidal 
damping of semimajor axis and random velocity namely 

v 2(vr) 2(vk)k 
atidai - ~r , {U) 

where r, v, and a are Cartesian position, velocity, and acceleration vectors, respectively (with 
r as the magnitude of the radial vector) and k is the unit vector in the vertical direction. 

In the code (to be discussed in §3), it is equations (9) through (12) which are used to 
incorporate the interaction between the embryos and planetesimals and the gas disc. 



3. N-body Methods 

We perform the numerical integrations using a parallel implementation of SyMBA (Dun- 
can et al. 1998) called miranda. SyMBA is a mixed-variable symplectic integrator based on 
the N-body map of Wisdom & Holman (1991) (see also Kinoshita et al. 1991) which has 
been improved to accurately handle close encounters between particles. As is done in the 
analytic model, we consider two classes of object, embryos (which gravitate and can accrete 
planetesimals) and field planetesimals (where we neglect all self-interactions.) We have made 
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several modifications to the original algorithm, mainly to include the effects of the gas disc. 
As introduced in §2.3, both embryos and planetesimals suffer the aerodynamic drag of eq. (9) 
and the embryos further suffer semimajor axis decay and random velocity damping accord- 
ing to the prescriptions of (10), (11), and (12). Accelerations resulting from the gas disc 
are treated as operators surrounding the (possibly recursively subdivided) drift operator in 
the basic democratic heliocentric step (Duncan et al. 1998) after the manner of Thommes 
(2001). However, in practice they are not applied during close encounters as the model for 
the gas density is almost certainly inapplicable in such situations, and we want to prevent 
artificial formation of tight binaries. Note that we neglect the gravitational potential due to 
the gas disc. 

One important approximation made in our N-body treatment of the field is the use of 
'super-planetesimals'. That is, in representing the field by a discrete number of particles, we 
use a large mass m sp for each object and imagine it represents an aggregate of underlying 
planetesimals of mass m. The aerodynamic drag that the object feels is that which the 
underlying m-mass object would, but when a merger occurs the mass m sp is used. Thus 
the integrated field objects are serving as dynamical tracers. Theoretical considerations 
(e.g. Binney & Tremaine 1987) show that the accretion rate of an embryo embedded in a 
planetesimal field should depend only the product of the number density of the planetesimals 
and their mass, provided that the collision timescale is much shorter than the embryo growth 
timescale. Numerical experiments varying the ratio m sp /m over orders of magnitude confirm 
that the accretion is insensitive to the exact value of m sp in our regime given that the 
characteristic embryo mass M ^> m sp and the number of super-planetesimals is large. 

It remains true that as a consequence of this approximation, at early times when the 
ratio of typical embryo mass M to m sp is at its lowest, we are inaccurately treating the 
effects of dynamical friction. However, the embryos are on circular, non-inclined orbits in 
any event due to the strong tidal eccentricity damping, and the aerodynamic drag is damping 
the planetesimals' random velocities, which suppresses the frictional stirring. As the drag 
becomes less effective due to the dissipation of the disc, the M/m sp ratio quickly increases, 
thereby reducing the inaccuracy. 

We adopt as starting conditions M/m sp > 5 (where M is the characteristic initial 
embryo mass). Note that this differs from the approximation often made by Kokubo and 
Ida where they instead increase the accretion radius of their particles by a factor /, thereby 
artificially speeding up their simulations by a factor f^ 13 for some (3 G [1, 2] (c.f. Kokubo & 
Ida 1996; e.g. Kokubo & Ida 2002). 

We also assume perfect accretion, which avoids the problem of collisional cascades (and 
the resulting increase in particle number, catastrophic in a nonstatistical code) at the cost 
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of inaccuracy in the mass spectrum; effects of this and other simplifications are discussed in 
§5. 

4. Simulations 

4.1. Generating initial conditions for the N-body simulations 

Since computer power is limited, we wish to begin our integrations as late in the second 
stage as possible. We will use the semi-analytic model of §2 to evolve the disc from our nomi- 
nal T=0 start until the number of embryos (i.e. the predicted number of embryos; recall that 
M(a) is continuous) is computationally tractable. This naturally produces a self-consistent 
gradient in embryo masses with the appropriate field population. We then discretize M(a) 
and £ m (a), yielding a set of embryos and a large number of super-planetesimal field parti- 
cles, typically ~8500 (subject to the constraints on mass discussed in the previous section), 
which constitute our N-body realization of the initial conditions for our simulations. 

To be specific, we set embryo and planetesimal densities pu = p m = 3.0 gem -3 , and 
take the underlying physical planetesimal mass to be m = 2.1-10" 6 M ffi (i.e. planetesimal 
radius r m = 100 km at 3.0 gcm~ 3 ). We set the efficiencies of aerodynamic drag and tidal 
eccentricity damping to unity, Co = l,c e = 1, and set rj = 0.001 at 1 AU. We take M , 
the initial embryo mass at T = for the semi-analytic model, to be 0.0015M e , which is 
much larger than m but two orders of magnitude below the likely final values M > 0.1M e . 
By introducing this larger M$ we advance the evolution (especially of the embryos at larger 
a) forward in time. This procedure is justifiable to the extent that the evolution we are 
bypassing is merely local accretion prior to any significant embryonic tidal migration. Thus, 
embryos in a simulation starting with smaller embryonic seeds would largely 'catch up' in 
mass to those in our simulations by the time significant migration sets in. 

We evolve this model until T = 0.3 Myr, which produces a satisfactory (~ 40 — 50) 
number of embryos, and is early enough that the use of a fixed b of 10 is acceptable. In this 
early phase, the runaway growth timescale (to return the system to oligarchy if the spacing 
is too large) and the scattering timescale (to do the same if the spacing is too small) are 
both short. Relaxation to the empirical equilibrium conditions begins almost immediately 
in the simulations, and we do not believe that our results are particularly sensitive to the 
fine details at the start. 

From this T = 0.3 Myr frame, objects are constructed from a = 0.75 AU to a x ks 2.5 
AU with embryo masses derived from M(a) (objects placed inwards, from 0.50 AU to 0.75 
AU, quickly fall off the inner edge of the simulation at 0.40 AU.) This sets the maximum 
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mass of a field object m sp (see §3) and from S m (a) the field objects can be built. Initial 
eccentricities and inclinations for embryos are set to 0.001 and 0.0005, respectively. Plan- 
etesimal eccentricities are drawn from a Rayleigh distribution of scale e eq using eq. (3), and 
we enforce % = e/2. (We observe that eq. (3) consistently predicts an equilibrium eccentricity 
which is slightly larger than the value to which the simulations rapidly relax. Experiments 
suggest that this is due to an overprediction of the stirring, whether the viscous or friction 
term, and not an underprediction of the effects of aerodynamic drag. This is consistent with 
the comments of Stewart & Ida (2000) on stirring efficiency.) All angles are randomized. 

As discussed above, the implications of the simplified models of §2 are that to preserve 
sufficient mass we will need to increase the surface density, and therefore (keeping in mind 
possible enhancements in the outer solar system) we focus on discs ~ 3 to 4 times the 
minimum mass disc. Specifications of the resulting conditions are summarized in table 4.1. 
We use values for Tdecay of 0.5, 1.0, and 2.0 Myr, and migration efficiency c a of 0.25, 0.50, 
and 1.00. For concreteness we choose the Tdecay = 1 Myr, c a = 0.5 case to discuss below. 

It was not obvious to what absolute time T the simulations would need to be advanced: 
after some experiments we chose a value of T = 20 Myr, when the gas has effectively vanished 
in all our simulations and the inter-embryo spacing is large enough that the timescale for 
mutual interaction is long. Note that we do not include Jupiter and Saturn in this simulation. 

For example, fig. 5 shows the initial (i.e. T = 0.3 Myr) conditions for the N-body 
simulation produced by the model for run C2. We see the clear decrease in embryo mass 
with increasing a characteristic of early times. In the innermost region, 0.70-0.80 AU, roughly 
50% of the material is contained in embryos, and in the outermost region (2.4-2.5 AU) this 
proportion drops to less than 5%. (We define / em b as the fraction of total mass in a region 
which is in the form of embryos.) 

The majority of the simulations were run at the Canadian Institute for Theoretical 
Astrophyics on the McKenzie parallel machine, a Beowulf cluster of ~ 256 dual-processor 
2.4 GHz Linux boxes, and the remainder were integrated on local machines at Queen's 
University. 
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Table 1: Simulation Parameters 



Name 


/cnh 


V 


Tdecay [Myr] 


C a 


N M 


N m 


CI 


3.0 


-1.0 


1.0 


0.25 


47 


7987 


C2 


3.0 


-1.0 


1.0 


0.50 


47 


7987 


C3 


3.0 


-1.0 


1.0 


1.00 


47 


7987 


C4 


3.0 


-1.0 


2.0 


0.25 


47 


7876 


C5 


3.0 


-1.0 


2.0 


0.50 


47 


7876 


C6 


3.0 


-1.0 


2.0 


1.00 


47 


7876 


C7 


3.0 


-1.0 
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Note. — / cn h and p arc defined through S so iid = (/enh ■ 7.1 gcm~ 2 )(a/AU) p where S so iid is the original 
surface density in solids at T=0. Td OC ay is the e-folding time of the gas dissipation, c a is the migration 
efficiency factor, and Nm and N m are the embryo and field particle counts. Each field particle is a super- 
planetesimal as described in §3. 
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4.2. Results 

As an example, time slices in (a, M) and (a, e) from C2, a representative run which 
demonstrated most of the typical behaviours, are displayed in figures 5 through 10. In this 
simulation, S m = 21.3 gem -2 (a/AU) _1 , c a = 0.5, and Td eC ay — 1 Myr. 

After 0.7 Myr, in figure 6, at T = 1 Myr (1 e-folding time of the gas decay), the system 
has undergone considerable evolution. Many embryos have undergone mutual mergers and 
some have migrated beyond the inner edge of the simulation and been removed. Figure 6 
shows wide ranges of roughly oligarchic behaviour where neighbouring objects have compa- 
rable mass. In the region from 0.5 to 1.6 AU, for example, most of the embryos have mass 
~ 0.18M e ± 0.02M e . At this time the region where f cmh ~ 0.5 extends from 1.0-1.5 AU. 
The variation of evolution timescale with semimajor axis is apparent, but both embryo-field 
and embryo-embryo mergers have occurred beyond 1.5 AU. 

In zones where an embryo has achieved some separation from its neighbours, 'Jacobi 
wings' can be formed (see Tanaka & Ida 1997). Jacobi wings are the characteristic config- 
urations - a curved V - produced in the (a, e) plane by an embryo embedded in a disc of 
planetesimals as it scatters them, so called because the shape of the scattering paths of the 
field particles are determined by the the conservation of the Jacobi integral in the restricted 
circular three-body problem. One such wing is apparent interior to the embryo near 0.7 AU 
in fig. 7, where the Jacobi wing has swept up a considerable amount of mass (0.2M®) during 
migration. 

At this time, 0.43M e of material has been removed from the simulation for approaching 
the Sun too closely (0.4 AU being the limit). Ultimately 2.1M® will be removed, and 90% 
of this is accomplished by 3 Myr. The probable fate of the removed material is discussed 
further in §4.5; there are good reasons to believe that most of the material which escapes 
makes it to the solar surface. With migration, in order to leave enough mass in the terrestrial 
region, it becomes necessary to make a substantial sacrifice to the Sun. 

At 2 Myr, in figure 7, when the gas has dropped to ~ 1/7.4 of its original density, more 
objects have been lost to the inner edge, and those that remain from 0.5 to 1.5 AU retain 
comparable mass (i.e. the deviation from the mean is less than a factor of two) and have 
grown to a mean mass just over 0.2M e . By 5 Myr (figure 8) the / emb < 0.5 region is now 
beyond 2.0 AU and we are approaching the end of the transition phase and the beginning 
of our end stage (where embryo-planetesimal accretion no longer plays a significant role.) 
The mean embryo mass in the 0.5-1.5 AU region is ~0.25M e , brought down somewhat by 
the objects from 1.3 to 1.4 AU. At this particular time the embryo spacing from 0.5 to 0.8 
AU is anomalously small, and this tightly-packed configuration quickly evolves to a wider 
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one. At this time, after 5 e-foldings, the gas density has dropped to ~ 1/150 of its original 
value and gas drag (whether aerodynamic or tidal) is not important. The corresponding 
increase in eccentricity pushes the accretion timescale up substantially, and the dominant 
growth mechanism will now be embryo-embryo mergers. 

Between 10 Myr and 20 Myr (figures 9 and 10) we arrive at a long-lived configuration 
of embryos: only one embryo merger occurs between 0.5 AU and 1.5 AU, and it occurs 
before 14 Myr. Evolution is continuing beyond 1.75 AU, and although limited, the remaining 
planetesimals are gradually being consumed. The mean embryo mass in the region of interest 
(0.5-1.5 AU) is ~0.45M e , and the spacing varies from b ~ 15 to b ~25, which is very well- 
separated. The eccentricities of the bodies (~ 0.02 — 0.04), which have been increasing over 
time, are now many times above their tidally-suppressed values. The region contains 2.58M e 
in the form of embryos and a residual planetesimal field of 0.23M e for a total of 2.81M e in 
the zone from 0.5-1.5 AU. 

As previously noted, roughly 2M e of material escapes the simulation to the inside. 
A substantial amount of mass remains beyond the zone of interest, however: 3.55M e of 
material is located beyond 1.5 AU. This is problematic, as there is nothing approaching a 
4M e object between the Earth and Mars; we return to this problem in §4.5. 

Viewing the embryo evolution in (M, a, t) space is instructive, as shown in fig. 11. Three 
different regimes are apparent: for a < 1.0 AU the embryos head directly for the Sun without 
scattering; for 1.0 AU < a < 1.5 AU the embryos 'slide' (i.e. collectively move in an orderly 
fashion, without orbit crossing) toward their final destinations, slowing as the gas decay 
lowers the migration rate; and for a > 1.5 AU the evolution is highly chaotic. Figure 12 
magnifies the innermost regions. 

At early times, when the embryo masses are low, the embryos can reorder themselves 
substantially and move several tenths of an AU without suffering a merger with another 
embryo. However, as time passes and the masses of the embryos grow (whether through 
accretion of planetesimals or direct embryo mergers), the migration pulls objects inwards 
at different rates. This leads to embryo segregration, and in some cases (typically where 
migration rates are high) to the formation of groups of two to four members which we call 
tidal convoys. (Although not mentioned explicitly, this effect also appears to have been 
present in the simulations of Papaloizou & Larwood 2000; see figure 6 in their paper.) 

At any time in these convoys the objects are typically ordered such that the mass 
increases outwards. This is self-selection: if an inner object is more massive than an outer 
one, differential migration will cause the objects to diverge, breaking the convoy. As one 
would expect, the migration rate of the convoys is above that of the less massive objects and 
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run C2 : time = 3.00e+05 [yr] 

o 




semimajor axis [AU] 



Fig. 5. — Initial conditions at T = 0.3 Myr for model C2. Large circles correspond to 
embryos (of area proportional to mass) and small circles to planetesimals; the lines through 
the large circles indicate a width of 10 Hill radii. The broken lines on the top graph indicate 
the amount of planetesimal material (in Earth masses) in a semimajor axis bin of width 0.1 
AU. 
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Fig. 7.— Configuration at T = 2 Myr for model C2. 
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run C2 : time = 5.00e+06 [yr] 



o 



CO 

d 
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Fig. 8.- 



Configuration at T = 5 Myr for model C2. 
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run C2 : time = 1.00e+07 [yr] 

o 



CO 

d 




semimajor axis [AU] 
Fig. 9.— Configuration at T = 10 Myr for model C2. 
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Fig. 10.— Configuration at T = 20 Myr for model C2. 
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Fig. 11. — Evolution of embryos with time in C2. The area of the circles scales linearly with 
embryo mass as indicated in the legend. For each embryo three lines are drawn, corresponding 
to the osculating values of the perihelion, semimajor axis, and aphelion. 
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run C2 : convoying and sliding regions 




time [Myr] 

Fig. 12. — The interior regions of C2 at early times. 
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below that of the most massive object as the outer objects 'push' on the inner ones. In many 
of our simulations, as in C2, these groupings are only apparent on the inner edges (where the 
object masses and gas density are highest) as the decay of the gas disc and corresponding 
weakening of the migration makes the convoys harder to form and maintain. If we set 
Tdccay = oo, tidal convoying becomes nearly inevitable, which raises serious questions about 
the applicability of analytic oligarchy models to the middle-to-late stage transition unless 
the gas is presumed to have vanished. 

In our simulations, while the gas remains, type I migration is remarkably efficient at 
capturing embryos into first-order mean motion resonances with each other, and the convoys 
usually show resonant relationships (two-object pairs are especially susceptible.) This en- 
hancement of resonance capture probability due to migration has been discussed in a different 
context by previous authors (see Peale & Lee 2002 for a possible explanation for the Laplace 
resonance among the Galilean satellites.) Most two- member tidal convoys are observed to 
be stable, but arguments to that end (e.g. Gladman 1993) are not easily generalized to the 
multiple-planet case, and in practice convoys involving three or more objects of different 
masses are eventually observed to merge. 

Beyond the convoying region, fig. 12 shows part of the 'sliding' region (and, conveniently, 
a pair showing intermediate behaviour.) As is apparent in fig. 11, when the masses and gas 
densities are such that the spacings will stay roughly constant during the migration, a group 
of embryos can move smoothly from further out in the disc suffering only a few mutual 
mergers. (Admittedly, in C2 the surface density of the gas disc varied as £ oc a -1 , and one 
can show (see section 4.5) that da/dt is constant for a fixed mass in this case. However, in 
the sliding region the masses vary by a factor of several, the gas density was decaying, and 
similar behaviours are seen for the £ oc cr 3 / 2 case. The effect is not a peculiarity of this 
particular density law, although it is a consequence of the near-constant migration rates that 
it and similar laws provide.) 

Finally, in the lowest-mass region, chaotic behaviour persists for the entire length of 
the simulation, although it does begin to settle down by the end. In one case, the embryo 
excursion reached 1.3 AU. Even in this chaotic regime, the broad outlines of the oligarchic 
model (roughly equally spaced embryos of comparable mass) were still respected, despite the 
fact that the simplest oligarchic picture is inapplicable as the embryos are being constantly 
reordered. This may imply that other mechanisms underly the oligarchic model's impressive 
robustness (Goldreich et al. 2004). 

These three regimes (convoying, sliding, and chaotic) were observed throughout the 
simulations - although their locations and strengths varied with mass, gas density, and 
stochastically-set spacing, as expected - and are likely frequent. 
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4.3. Final Characteristics 

Time slices of the final configurations for simulations C1-C14 are plotted in figures 13, 
14, and 15. 

The resulting mass in the 0.5-1.5 AU region in the form of embryos ranges from 0.81M e 
to 3.27M ffi , and for planetesimals from 0.14M ffi to 0.44M ffi . Some of the runs (C2, C3, C4, 
C6) yield very oligarchic-looking outcomes where the variation in mass between embryos 
is low, and the characteristic mass varies from 0.2 to 0.4 M e . Others (CI, C7, C8, C9, 
CIO) produce systems with embryos of mass approaching lM e , at least some of which are 
candidate planets. The remaining simulations are difficult to classify, although most show 
regions oligarchic in appearance with a few outliers (e.g. C12, C14.) C5 produced the least 
mass (embryo mass 0.81M©, field mass 0.44M ffi ) for reasons to be discussed in §4.6, and CI 
the most (3.27M e , 0.16M®). By accident, CI produced two dominant objects of terrestrial 
mass - Venus and Earth analogues. 

Figures 16 and 17 show the resulting amounts of embryo and field mass for all simulations 
at 1, 10, and 20 Myr. Most simulations show an increase in embryo mass with time, and 
a substantial number of simulations landed in our nominal target region of 1.5 to 2.5M e in 
embryo mass and < 0.5M ffi in planetesimal mass. 

The original expectation was for all the simulations to produce oligarchic-looking results 
(as in C3), not direct planet production (as in CI), but it is not surprising that with more 
massive embryos at least some embryo-embryo mergers would occur as the gas vanished, 
and it takes only a few such mergers at late times to produce a substantial planet. Of 
the fourteen simulations, four produced less than two Earth masses (our nominal target) of 
material in the region and the remaining ten produced more, leading one to suspect that 
our enhancements may have been too high. However, this depends entirely on how much 
of the material is removed during whatever processes turn the resulting configurations into 
terrestrial systems: they are quite widely spaced. 

4.4. Evolution of inter-embryo spacing 

As noted by Goldreich et al. (2004), the argument offered in Kokubo & Ida (1998) to 
derive an estimate for the inter-embryo spacing b as a function of M and £ is of dubious 
applicability. Recall that the spacing b in their model is derived by equating the decrease in b 
due to embryo growth (due to planetesimal accretion) with an increase in b due to scattering. 
The two-body scattering formula of Petit & Henon (1986) has a problem in the oligarchic 
context: namely, there are embryos on both sides of our scattering pair who are instead 
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Fig. 13. — Final configurations in mass and semimajor axis for runs C1-C6. M T and m T 
indicate (in Earth masses) the total mass in embryos and in planetesimals, respectively, in 
the region [0.5 AU, 1.5 AU]. 
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Fig. 14.— Same as fig. 13, but for runs C7-C10. 
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Fig. 15. — Same as fig.13, but for runs C11-C14. 
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Fig. 16. — Embryo mass versus planetesimal mass within the 0.50-1.50 region at 1, 10, and 
20 Myr for runs C1-C7. The nominal target zone (from 1.5 to 2.5 M e in embryo mass and 
< 0.5M e in planetesimal mass) is shaded. 
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Fig. 17.— Same as fig.16, but for runs C8-C14. 
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trying to scatter the two objects closer together. Naively one would expect this to turn the b 
growth from a monotonic process into a much slower random walk. The formula also assumes 
that the eccentricities of the two objects are near zero at conjunction, an assumption which 
breaks down in our model as the damping timescale increases with the dissipation of the gas. 
(Simple experiments confirm both suspicions: the multiple-embryo case shows substantially 
slower b growth than the two-embryo case, and introducing a decay in the strength of the 
eccentricity damping in the two-embryo case also washes out the growth.) Finally, when 
type I migration is active, differential migration of equal-mass objects will tend to increase 
b, and this is a non-negligible effect in our regime. 

Nevertheless, our numerical experiments suggest that b ~ 10 is a reasonable approxi- 
mation during the period we apply the semi-analytic model, possibly accidentally: embryos 
quickly scatter when b < 5 and growth is a weak function of M (b oc M 2//15 in Kokubo 
& Ida 1998) so substantial masses are required before b = 10 becomes unacceptable. In 
any case, corrections to the 'wrong' chosen b occur on short timescales, and during testing 
no significant correlations between initial b (granted b ~ 10) and the results were observed, 
suggesting that the equilibration process leaves little signature. 

Figure 18 shows the evolution of the mean mass-weighted b in the region a = 0.50 — 
1.50AU for all simulations. Although the specific evolution varies from integration to in- 
tegration, the trend of increase in b from ~ 10 to ~ 20 is clear. The trend to large b is 
accentuated by the tidal migration: in our control runs without tidal migration (not shown) 
we find values of b ~ 10 — 15 were common in the late stages. The spikes near 4 and 6 Myr 
are the result of unusual structures being formed; the 4 Myr spike in run C5 which reaches 
off the scale is to be further discussed in §4.6. One concern is that we find spacings of more 
than 20 Hill radii are typically sufficient to produce stable systems (at least on timescales of 
100 Myr, the expected formation timescale). Preliminary simulations suggest that the intro- 
duction of Jupiter and Saturn stirs the system enough to produce embryo interaction, and 
resonance sweeping involving the removal of the disc potential (here neglected: see Nagasawa 
et al. 2003) is another possibility. We will return to this issue in a subsequent paper. 



4.5. Material Transport 

The gas drag (both aerodynamic and tidal) tends to drive embryo and field material 
inwards. Scattering amongst embryos, or between embryos and field objects, can counteract 
this and move material outwards. The question of how the transport processes affect the 
resulting proportions of mass in the final embryos therefore arises (i.e. what amount of 
mass, from where.) Figures 19 and 20 show the resulting material fraction by source region 
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Fig. 18. — Evolution of mean inter-embryo spacing measured in units of single-planet Hill 
radii in region 0.5-1.5 AU. 
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( whether embryo or planetesimal.) 

Substantial amounts of the material which ends up in the 0.5-1.5 AU region comes from 
beyond 1.5 AU, and even from beyond 2 AU. In some rare cases (e.g. the embryo at 1.2 
AU in CI, or at 0.87 AU in C8), an embryo in fact consists mostly of material from beyond 
2 AU. It is not surprising there are general trends with rd cca y Indeed, in fig. 19, we see 
in runs C4, C5, and C6 (each with Td eC ay — 2 Myr), more material from beyond 2 AU is 
incorporated than in C1-C3 (with decay timescale 1 Myr.) A somewhat weaker trend with 
migration efficiency c a is also evident. Comparatively little material originally from 1 AU 
becomes part of embryos beyond 1.2 AU. We do not believe that our lack of embryos and 
planetesimals beyond 2.5 AU at T = 0.3 Myr plays a significant role in these results: in the 
majority of the resulting embryos, material from beyond 2.25 AU contributes little, and the 
arguable exception of C6 is the simulation with the most migration (the result of a 2 Myr 
dissipation timescale and an efficiency c a = 1). 

The large degree of inward mixing suggests that present-epoch local chemistry may 
be probing much further out in the planetary disc than might have been expected without 
migration, at least if countervailing effects such as disc turbulence are limited. Consequences 
of different models of material transport during terrestrial planet formation on presently 
observed water abundances are discussed in Lunine et al. (2003) and Raymond et al. (2004). 

Recently Weidenschilling, in unpublished work, has also been investigating terrestrial 
formation in enhanced discs with migration using a hybrid dynamical and statistical tech- 
nique and reached broadly consistent conclusions (private communication, DPS 2004.) How- 
ever, he reports a definite tendency for material to accumulate in the inner regions. 

As we remove objects when they get below 0.4 AU, and do not simulate any objects 
inside of 0.75 AU at our start of 0.3 Myr, we are insensitive to the production of a terrestrial- 
mass Mercury. This decision was motivated by our early experiments, which demonstrated 
that the oligarchic predictions of the §2 model were reliable (except for certain details of the 
b spacing) over a large mass range before type I migration played a role and were tolerable, 
although less accurate, at tidally significant masses. The resulting embryos were massive 
enough that they invariably migrated inside of 0.4 AU before the gas vanished. Although 
these early embryos escaped our region of interest, this is not a guarantee that they would 
succeed in migrating all the way to the Sun. To investigate the likely future of embryos 
whose dynamics we stopped tracing, we can take the objects removed for falling off the inner 
edge of our simulation and integrate the tidal migration equations (8) and (7). For £ oc a -1 , 
then da/dt oc exp(— t/Td eca y)) and for X oc a~ 3//2 , da/dt oc a -1 / 2 exp(— t/rd eca y)- 

Figure 21 shows the projected evolution of the escaped embryos, considered one at a 
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time. The majority of objects (80%) succeed in making it to the solar radius (~ 0.005 
AU), and so our removal of them from the simulation is defensible. In over half of the 
simulations, no embryo which was removed would have survived. In several runs (CI, C2, 
C4, C9, C10, Cll), multiple embryos would have survived, although in only four of these 
(C2, C4, C10, Cll) would an embryo have survived inside of 0.2 AU. Note this analysis is 
restricted to the embryos which were actually integrated and then removed, not to embryos 
which should have been there but never were (e.g. embryos initially located at 0.5 AU). A 
fortiori most of the embryos we did not include inside of 0.75 AU would have also made it 
to the surface of the Sun, at least for the 1 Myr and 2 Myr decay timescales. In the 0.5 
Myr case, with its resulting lower migration, although we are not missing much material 
which should end up in the 0.5-1.5 AU region, we could be missing a substantial amount 
between 0.2 and 0.4 AU by our choice of initial embryo range. As for planetesimals, given 
the effectiveness of migrating Jacobi wings as a sweeping mechanism, it seems likely that 
many of the planetesimals not incorporated into the interior embryos would move with the 
flow, and given the strong dependence on semimajor axis of formation timescale there are 
unlikely to be many left at 0.4 AU when tidally-significant embryos are emerging at 1 AU 
and beyond. In summary, for r dccay > 1 Myr, we expect that we are missing only a few 
embryos, and do not believe our choice of inner edge significantly affects our conclusions. 

The outer regions present their own set of difficulties: an average of 3M e of material is 
left beyond 1.5 AU (recall that the initial outermost embryo in each simulation is located 
at ~ 2.5 AU), which is an order of magnitude above the mass of Mars (0.11M ffi ). As is 
evident from the discussion of §4.3, we have provided an upper bound on the necessary disc 
enhancement to survive type I migration, as our enhancements consistently produce more 
mass than desired. Decreasing the enhancement to reduce this overproduction will mitigate 
this problem somewhat, but even a one-third decrease will still leave ~ 2M e of material. 

This problem is not unique to our migration models, however: the naive minimum mass 
model, Sgoiid = 7.1 gcm~ 2 (r/AU)~ L5 , itself produces ~1.1M® of material between 1.5 and 
2.5 AU. The majority of the unwanted material in this region in our simulations is simply a 
reflection of how much was there to begin with. One can deal with this problem in several 
ways: assume a primordial decrease in original surface density in the region between Earth 
and Mars; increase the material transport rates in some fashion, e.g. collisional grinding 
producing large amounts of dust which leaves the region very quickly due to aerodynamic 
drag; and so on. Any approach which solves the problem for non-enhanced models without 
migration will have a natural analogue in our enhanced models with migration, and therefore 
the problem is no worse in our model than for the standard scenario. 

Of course, it is unlikely that a simple power-law profile accurately describes the density 
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Fig. 21. — Projected future evolution of embryos removed from the simulation 
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of the protoplanetary disc (either solids or gas) to within a few radii of the Sun. If the solid 
density reaches a local maximum near 1 AU, it may be possible to avoid any interior and 
exterior problems entirely. 

4.6. Rings 

An unexpected side effect of neglecting the self-interactions of the field particles became 
apparent during some of the simulations. Occasionally, due to an early embryo-embryo 
merger, an object is produced which is prematurely massive relative to its neighbours. The 
increase in mass results in a higher migration rate than the next outer embryo, and a gap 
opens up between them where the only objects are planetesimals. This region devoid of 
embryos can also be created by a major scattering event. In a weak sense this process occurs 
during every merger, but most of the time the gap is quickly filled by another embryo, 
whether through migration or scattering. However, under certain circumstances, this gap is 
not closed quickly with the entrance of an embryo from further out in the disc, but instead 
persists. In the absence of embryos there is nothing to stir the field, and the aerodynamic 
drag decreases the field eccentricity, resulting in a thin ring of planetesimals (except where 
an exterior embryo arrives to form a new Jacobi wing.) In extreme cases the width of the 
ring is substantial (e.g. 0.2 AU), corresponding to a large amount of mass (~1.0M e ). At this 
point it can serve as a wall, and act as a barrier to the entry of inward- migrating embryos 
from beyond the ring. 

There are several possibilities for the resulting behaviour. Sometimes the ring is quickly 
disrupted, especially if it is small. At other times the ring endures until an exterior embryo 
pushes through. This often leads to a surprisingly quick transition of the embryo through 
the ring, a consequence of the process of forming Jacobi wings. The ring may or may not 
survive this embryo crossing. Sometimes the ring endures, and the evolution stalls as the 
embryo masses reached exterior to the ring are insufficient to push through the barrier. 

An example of an extreme scenario (from run C5) is demonstrated in figure 22, the 
most massive (by a factor of several) and longest-lived (also by a factor of several) ring in 
our simulations. At 2 Myr, the system is unremarkable: between 1.0 and 1.7 AU there are 
four embryos separated by ~0.1 AU, with planetesimals throughout. By 2.5 Myr, however, 
scattering events have resulted in a region from 1.3 AU to 1.6 AU without any embryos, but 
with almost an Earth mass of planetesimal material. After another half million years, the 
innermost embryos have migrated out of the region, and only one embryo of ~0.2M e remains 
to control the dynamics. (This wide embryo-free region is responsible for the highest spike 
in b in fig. 18.) Over the next 8 Myr, the ring moves from ~1.5 AU to ~1.25 AU as a result 
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of both aerodynamic drag (over the first several million years, at any rate; the dissipation 
timescale in C5 is 2 Myr) and the more significant push given by the exterior embryo. During 
this push, the embryo consumes enough of the ring during the Jacobi scattering to grow from 
0.2 to 0.3M e . Finally, an encounter between the shepherd embryo and and its neighbour 
scatters the shepherd into the ring by 11.5 Myr. When the embryo enters, it quickly migrates 
through and disrupts the ring. By 12 Myr the ring is gone. 

We emphasize that this is an extreme case, and that when the rings do not form or are 
quickly disrupted by the embryos (before field interactions would be expected to do so) the 
two-component embryo/field approximation we make remains acceptable. 

Could these rings be physical, or are they merely artifacts of our simplified treatment 
of the planetesimal field? It is true that some early mergers and scattering events are 
going to occur, and differential migration will then attempt to produce a gap. However, 
in a system with a more realistic mass spectrum (with objects of intermediate mass) the 
second-tier objects which are not accreted during the migration of the overmassive embryo 
might ascend to be the new oligarchs in the absence of competition. Furthermore, eventually 
mutual planetesimal-planetesimal encounters (which we neglect) should produce new seed 
objects for oligarchic growth within the ring. Rings which last for less than the timescale 
for internal disruption seem possible, and may occur in real planetary systems. In this case, 
their chief effect would seem to be slowing the migration process for embryos with exterior 
semimajor axes. 



5. Discussion 

Various authors, looking with dismay at the unhappy consequences of type I migration 
for forming giant cores, have concluded that the effect simply does not occur. Some recent 
work (Menou & Goodman 2004) has also argued that the timescale for migration in some 
circumstances may be substantially larger than we have assumed here. This may ultimately 
prove correct, but several points are worth making in response. First, given the complexity 
of the local gas physics near the embryo, it is probably premature to rule out near-nominal 
migration; a clear understanding of the behaviour will need to wait for the next generation 
of hydrodynamic simulations. For example, it is entirely possible that whatever effects 
hypothetically shut off type I migration before the opening of a gap only become significant 
for embryo masses near or beyond Earth mass. If so, the intuition that type I is migration 
is too destructive to have occurred in our system could be completely justified regarding 
the giant planets but not relevant for terrestrial formation. Disc properties such as opacity 
on which migration may sensitively depend are extremely poorly constrained, and saying 
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Fig. 22. — Extreme example of the formation of a planetesimal ring. The ring endures from 
2.5 Myr to 11.0 Myr in run C5. The left axis shows the eccentricity of the objects (dots) 
while the right axis shows the mass of planetesimals in bins of width 0.1 AU (solid lines). 
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anything firm will require improved proplyd observations. 

That said, it seems possible that the gas disc is far less smooth than we have assumed, 
and type I migration may be halted by abrupt changes in gradient, opacity, or density at 
unknown locations. For example, Matsuyama et al. (2003) develops models for the pho- 
toevaporation of the disc which result in a zone of depleted gas which could protect the 
terrestrial region from marauding proto- Jovian embryos. Although these specific models 
have dynamical consequences which make them difficult to accept at face value, the possi- 
bility of significant disc inhomogeneities causing similar effects remains. Unfortunately all of 
this will clearly be very difficult to model. Disc turbulence can wash out the effect of tidal 
migration on short timescales (Laughlin et al. 2004; Nelson & Papaloizou 2004) but may not 
remove a net drift; we hope to study oligarchy in a turbulent disc in future work. 

Several compromises made to complete the integrations with the available resources may 
affect the results. The neglect of field self-interaction contributes to the formation of rings, 
although to the degree that the behaviour of the system truly is oligarchic, the embryo- 
field model should be a reasonable approximation. Nevertheless, our representation of the 
mass spectrum is limited, and no field object can promote itself to an embryo. We have 
chosen 100 km as the underlying planetesimal size, which is both plausible (the formation 
timescale for such objects is comparable to the absolute time when we begin the simulation) 
and conventional in terrestrial formation simulations. Smaller planetesimal particles, such as 
those produced in a fragmentation spectrum, would lead to different accretion rates (Rafikov 
2004) in potentially complicated ways, and would also be likely to result in higher migration 
rates for planetesimals due to aerodynamic drag, possibly requiring greater disc enhancement. 

The artificial increase of the seed mass for the embryos may also cause problems. It 
changes the gradient in embryo mass, and means that at a given time, the outermost embryos 
are more massive relative to their inner neighbours than they should be. When the evolution 
is local, this is defensible, but during migration it will lead to an overestimation of the 
transport. The overestimation should be mild in practice as the outermost objects are sub- 
migratory at T = 0.3 Myr, when our N-body simulations begin. 

6. Conclusion 

We have investigated the likely conditions at the beginning of the terrestrial endgame 
as a result of type I migration during the mid-to-late transition. For reasonable values of the 
parameters we find that plausible progenitors for the terrestrial planets can survive despite 
type I migration at near the nominal rate, provided the disc density is enhanced above the 
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minimum model. In particular, in order to get the right amount of material left in the 
terrestrial region with a gas dissipation timescale of ~ 1 Myr, one must work with discs 
several times more massive than the minimum model, although a lower enhancement than 
we used here would likely suffice. (Preliminary investigations to be reported elsewhere using 
discs of ~ 2 times the minimum mass model are encouraging, and show good agreement 
with expectations.) As a consequence, the resulting systems at the beginning of the late 
phase show significant differences from the conditions normally considered (e.g. Chambers 
& Wetherill 1998; Chambers 2001): (1) there is often a large amount, ~0.40M e of material 
remaining in the form of planetesimals (supporting the direction taken in Chambers 2001) 
at the time of dissipation; (2) the embryos tend to have larger masses, with mode 0.4M e ; 
and (3) the embryos are separated by ~ 20 — 25 single-planet Hill radii. 

In a forthcoming paper we take the next step and study what must be done to turn our 
new endgame conditions into terrestrial planets. As our final configurations are stable when 
left alone, we will require the introduction of new events (e.g. the formation of Jupiter and 
Saturn as in Kominami & Ida 2004) - to complete the terrestrial construction project. 
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